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Hydrogen and deuterium chemisorption on a single layer of graphene has been studied by path- 
integral molecular dynamics simulations. Finite-temperature properties of these point defects were 
analyzed in the range from 200 to 1500 K, by using a tight-binding potential fitted to density- 
functional calculations. On one side, vibrational properties of the adatoms are studied at their 
equilibrium positions, linked to C atoms. The vibrations display an appreciable anharmonicity, 
as derived from comparison between kinetic and potential energy, as well as between vibrational 
energy for hydrogen and deuterium. On the other side, adatom motion has been studied by quantum 
transition-state theory. At room temperature, quantum effects are found to enhance the hydrogen 
diffusivity on the graphene sheet by a factor of 20. 
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I. INTRODUCTION 

In recent years there has been a surge of interest in 
carbon-based materials. Among them, those formed by 
C atoms with sp 2 hybridization have been intensively 
studied, as is the case of carbon nanotubes, fullerenes, 
and graphene. The latter, in particular, is up to now the 
only real two-dimensional crystal, with exotic electronic 
properties X& 

Carbon-based systems, in general, are considered as 
candidates for hydrogen storage 4 Also, chemisorption on 
two-dimensional systems, such as graphene, can be im- 
portant for catalytic processes.— The interest on hydrogen 
as an impurity in solids and on surfaces is not new, and 
dates back to many years. Even though this is one of 
the simplest impurities, a thorough understanding of its 
physical properties is complex due to its low mass, and 
requires the combination of advanced experimental and 
theoretical methods^^ 

Experimental investigations on atomic, isolated hydro- 
gen on graphene have been so far scarce, since H is diffi- 
cult to detect. Recently, it has been clearly observed by 
transmission electron microscopy, and its dynamics were 
analyzed in real timeX In general, apart from its basic 
interest as an isolated impurity, an important property of 
hydrogen in solids and surfaces is its ability to form com- 
plexes and passivate defects, which has been extensively 
studied in the last twenty years^^ 

From a theoretical viewpoint, atomic hydrogen on 
graphene has been studied by several authors using 
ab-initio methods i 4 ^ ' 10 ' 11 ' 12 ' 13 It is generally accepted 
that chemisorption of a single hydrogen atom leads to 
the appearance of a defect-induced magnetic moment 
on the graphene sheet, along with a large structural 
distortion^ ' 10 ' 11 However, standard electronic-structure 
calculations, in spite of their quantum mechanical char- 
acter, usually treat atomic nuclei as classical particles, 
and typical quantum effects as zero-point vibrations are 
not directly available from the calculations. Such quan- 
tum effects may be relevant for vibrational and electronic 



properties of light impurities like hydrogen, especially at 
low temperatures. 

Finite-temperature properties of hydrogen-related de- 
fects in solids have been studied by ab-initio and tight- 
binding (TB) molecular dynamics simulations. In many 
earlier applications of these methods, atomic nuclei were 
treated as classical particles! 14 ' 15 : 16 To consider the quan- 
tum character of the nuclei, the path-integral molecular 
dynamics approach results to be particularly suitable. 
In this procedure all nuclear degrees of freedom can be 
quantized in an efficient way, allowing one to include 
both quantum and thermal fluctuations in many-body 
systems at finite temperatures. Thus, molecular dynam- 
ics sampling applied to evaluate finite-temperature path 
integrals allows one to carry out quantitative studies of 
anharmonic effects in condensed matteri 17 ' 18 

In this paper, the path-integral molecular dynamics 
(PIMD) method is used to investigate the role of the im- 
purity mass on the properties of hydrogenic point defects. 
We study isolated hydrogen and deuterium (D) on a 
graphene sheet. Special attention has been laid upon the 
vibrational properties of these impurities, by consider- 
ing anharmonic effects on their quantum dynamics. The 
results of these calculations show that such anharmonic 
effects lead to an appreciable deviation of the vibrational 
energy of the impurities, as compared to a harmonic ap- 
proximation. Also, the atomic diffusion is found to be en- 
hanced with respect to the classical limit, for both H and 
D. Path- integral methods analogous to that employed in 
this work have been applied earlier to study hydrogen 
in metal s 17 : 19 ' 20 and semiconductors. 21 ' 22 : 23 ' 24 In connec- 
tion with the present work, hydrogen has been stud- 
ied inside and on carbon nanotubes by diffusion Monte 
Carlo^ 5 -^ 

The paper is organized as follows. In Sec. II, we de- 
scribe the computational methods employed in our cal- 
culations. Our results are presented in Sec. Ill, dealing 
with the energy of the defects, vibrational properties, and 
impurity diffusion. In Sec. IV we summarize the main re- 
sults. 
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II. COMPUTATIONAL METHODS 

In this Section we present the computational methods 
employed in our simulations. On one side, in Sec. II. A 
we introduce the PIMD method used to obtain equilib- 
rium properties related to the hydrogenic defects. On the 
other side, in Sec. II B we discuss a procedure to calcu- 
late rate constants for impurity jumps, in the context of 
transition-state theory. 



A. Path-integral molecular dynamics 

Our calculations are based on the path-integral formu- 
lation of statistical mechanics, which is a powerful non- 
perturbative approach to study many-body quantum sys- 
tems at finite temperatures. In this approach, the parti- 
tion function is evaluated through a discretization of the 
density matrix along cyclic paths, composed of a num- 
ber L (Trotter number) of "imaginary-time" stepsj 27 ' 28 
In the numerical simulations, this discretization gives rise 
to the appearance of L "beads" for each quantum par- 
ticle. Then, this method exploits the fact that the par- 
tition function of a quantum system can be written in 
a way formally equivalent to that of a classical one, ob- 
tained by substituting each quantum particle by a ring 
polymer consisting of L classical particles, connected by 
harmonic springs i 17 ' 18 Here we employ the molecular dy- 
namics technique to sample the configuration space of 
the classical isomorph of our quantum system (N car- 
bon atoms plus one impurity). Calculations were car- 
ried out in the canonical ensemble, using an originally 
developed software, which enables efficient PIMD simu- 
lations on parallel computers. The algorithms employed 
to integrate the equations of motion were based on those 
described in the literature ! 29 ' 30 

The calculations have been performed within the adia- 
batic (Born-Oppenheimer) approximation, which allows 
us to define a potential energy surface for the nuclear co- 
ordinates. An important question in the PIMD method 
is an adequate description of the interatomic interac- 
tions, which should be as realistic as possible. Since 
employing true density functional or Hartree-Fock based 
self-consistent potentials would restrict enormously the 
size of our simulation cell, we have derived the Born- 
Oppenheimer surface for the nuclear dynamics from 
an efficient tight-binding effective Hamiltonian, based 
on density functional (DF) calculations^ The capabil- 
ity of TB methods to reproduce different properties of 
molecules and solids was reviewed by Goringe et al£^ In 
particular, the reliability of our TB Hamiltonian to de- 
scribe hydrogen-carbon interactions in carbon-based ma- 
terials has been checked in previous work ; 24 ' 33 and ac- 
cording to those results we expect that systematic errors 
in calculated diffusion barriers for H in these materials 
are less than 0.1 eV. An advantage of the combination 
of path integrals with electronic structure methods is 
that both electrons and atomic nuclei are treated quan- 



tum mechanically, so that phonon-phonon and electron- 
phonon interactions are directly taken into account in the 
simulation. 

Simulations were carried out in the NVT ensemble on 
a 4 x 4 graphene supercell of size 4a = 9.84 A with pe- 
riodic boundary conditions, containing 32 C atoms and 
one adatom. For comparison, we also carried out sim- 
ulations of graphene without impurities, using the same 
supercell size. We have checked that larger supercells, 
i.e., 5x5 give within error bars the same results as those 
derived below from the 4x4 supercell. Sampling of the 
configuration space has been carried out at temperatures 
between 200 and 1500 K. For a given temperature, a 
typical run consisted of 2 x 10 4 PIMD steps for system 
equilibration, followed by 10 6 steps for the calculation of 
ensemble average properties. To have a nearly constant 
precision in the results at different temperatures, we took 
a Trotter number that scales as the inverse temperature, 
so that LT = 6000 K. For comparison with the results of 
our PIMD simulations, we have carried out some classi- 
cal molecular dynamics simulations with the same inter- 
atomic interaction, which is achieved by setting L = 1. 
The quantum simulations were performed using a staging 
transformation for the bead coordinates. Chains of four 
Nose-Hoover thermostats were coupled to each degree of 
freedom to generate the canonical ensemble ^ The equa- 
tions of motion were integrated by using the reversible 
reference system propagator algorithm (RESPA), which 
allows us to define different time steps for the integra- 
tion of the fast and slow degrees of freedom^ The time 
step At associated to the DF-TB forces was taken in the 
range between 0.2 and 0.5 fs, which was found to be ap- 
propriate for the atomic masses and temperatures studied 
here. For the evolution of the fast dynamical variables, 
that include the thermostats and harmonic bead interac- 
tions, we used a time step St = At/4. More details on 
the actual implementation of the simulation method can 
be found elsewhere i 3 - 3 -^ 



B. Quantum transition-state theory 

Classical transition-state theory (TST) is a well- 
established method for calculating rate constants of in- 
frequent events. An important element in this compu- 
tational method is the ratio between the probability of 
finding the system at a barrier (saddle-point) and at its 
stable configuration. A quantum extension of this theory 
has been developed in the context of path integrals, with 
the purpose of studying the kinetics of processes involving 
light atoms i 17 ' 36 ' 37 This quantum approach allows one to 
relate the jump rate k to the probability density of the 
center-of-gravity (centroid) of the quantum paths of the 
jumping atom, defined as 

i=l 
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Xj being the coordinates of the "beads" in the associated 
ring polymer. Then, k is related to the ratio P c between 
the equilibrium probability of finding the centroid at a 
saddle-point (say x*) and at a stable site (say Xo)i 36 ' 38 
Namely: 



where v is a factor weakly dependent on temperature, 
taken to be the thermal velocity v = ^2/(n(3m) of the 
jumping impurity, and / is the distance between Xo and 
x*. Note that, apart of typical quantum effects, the 
1/y/m factor in v will favor a faster jump rate of the 
lighter atoms, as in classical TST. The probability ratio 
P c can be written as exp(— 0AF), AF being an effective 
free-energy barrier, given by the reversible work done on 
the system when the impurity centroid x moves along a 
path from xo to x*: 

AF = - ( f(x)dx, (3) 

where f (x) is the mean force acting on the impurity with 
its centroid fixed on x at temperature T: 

f(x) = -<V x F(R))x. (4) 

Here V(R) is the potential energy, R being in our case a 
3(A^ + l)-dimensional vector (N carbon atoms plus one 
impurity). The average value in Eq. (j4|) is taken over 
quantum paths with the centroid of the impurity fixed on 
x at temperature T. Thus, the jump rate (a dynamical 
quantity) is related to the free energy difference AF (a 
time-independent quantity), so that its temperature de- 
pendence can be obtained from equilibrium simulations 
without any direct dynamical information* 3 ^ Quantum 
effects that may give rise to substantial deviations from 
the classical jump rate are taken into account within this 
kind of quantum TST. In this way, jump rates for kinetic 
processes can be obtained for realistic, highly nonlinear 
many-body problems. The reliability of this method to 
calculate free-energy barriers and jump rates has been 
discussed in Refs. Il9ll38l . In particular, it was argued 
that this method can be inaccurate in the presence of 
asymmetric barriers at low temperatures. In our case, 
the error bar of the calculated barriers is expected to be 
less than the intrinsic error of the employed tight-binding 
potential. 

In this kind of simulations, adatom motion was re- 
stricted along a reaction coordinate defined by the cen- 
troid of the quantum path. Note, for comparison, that in 
the equilibrium PIMD simulations described in the pre- 
vious section, there is no restriction on the motion of the 
H (D) centroid. The force f(x) has been evaluated at 
11 points along the reaction coordinate connecting the 
lowest-energy configuration (H linked to a C atom) with 
the saddle point of the energy surface, and then the inte- 
gral in Eq. ([3]) was calculated numerically. For each point 
in the integration path, we generated 5000 configurations 



for system equilibration, and 2 x 10 4 configurations for 
calculating the mean force at a given temperature. More 
technical details can be found elsewhere . 17 ' 40 ' 41 



III. RESULTS 

A. Vibrational energy 

We first discuss the lowest-energy configuration for the 
hydrogenic impurities on a graphene sheet, as derived 
from classical calculations at T = 0, i.e., point atomic 
nuclei without spatial derealization. The impurity binds 
to an C atom, which relaxes out of the sheet plane by 
0.46 A, with a bond distance between C and impurity of 

1.17 A. These results are in line with those reported in 
the literature, and in particular with the breaking of a 
7r bond and producing an additional a bond, changing 
the hybridization of the involved C atom from sp 2 to 
sp 3 ,! ' 10 ' 11 Assuming the host C atoms fixed in the relaxed 
geometry, one can calculate vibrational frequencies for 
the impurity in a harmonic approximation. Thus, we find 
for hydrogen a frequency u>± = 2555 cm -1 for stretching 
of the C-H bond (perpendicular to the graphene sheet), 
and uj\\ = 1186 cm -1 for vibrations parallel to the plane 
(twofold degenerate). 

We now turn to the results of our simulations at fi- 
nite temperatures, and will discuss the energy of the 
hydrogenic defects.. The internal energy of the system 
graphene plus impurity, E(T), at temperature T can be 
written as: 

E(T) = E min + E V {T) , (5) 

where E m i n is the potential energy for the classical mate- 
rial at T = (point-like atoms on their equilibrium posi- 
tions), and E V (T) is the vibrational energy of the whole 
system. Then, E V (T) can be obtained by subtracting the 
energy E m { n from the internal energy derived from PIMD 
simulations. In Fig. 1 we show the temperature depen- 
dence of the vibrational energy E v for a 4 x 4 graphene 
supercell including an H atom (circles). For comparison 
we also show E v for a pure graphene sheet (squares) . At 
300 K, the vibrational energy of graphene amounts to 

6.18 eV per simulation cell, i.e., 0.19 eV/atom. As ex- 
pected, E v increases as temperature is raised, and even- 
tually converges to the classical limit E^ 1 — ZNksT at 
high T. 

An interesting characteristic of the different hydrogenic 
defects (H or D) is their associated vibrational energy. At 
a given temperature, this energy is defined as the differ- 
ence AE V = E V (32G + Imp) - E V (32C), where 'Imp' 
stands for H or D. Note that AE V defined in this way 
is not just the vibrational energy of a given adatom on 
graphene, but it also includes changes in the vibrations 
of the nearby C atoms. Shown in Fig. 2 is the vibra- 
tional energy associated to the hydrogenic adatoms as 
a function of temperature. Symbols indicate results of 
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FIG. 1: Temperature dependence of the vibrational energy 
E v of the 4x4 graphene supercell with one hydrogen, as 
derived from PIMD simulations (circles). For comparison, we 
also present results for a pure graphene supercell (squares). 
Dashed lines are guides to the eye. Error bars are in the order 
of the symbol size. 
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FIG. 2: Vibrational energy AE V of the hydrogenic defects 
as a function of temperature. Results are shown for hydro- 
gen (squares) and deuterium (circles). For comparison, we 
also present the vibrational energy obtained in a harmonic 
approximation with frequencies wy and u>± corresponding to 
H and D (solid lines). Dashed lines are guides to the eye. 



PIMD simulations for H (squares) and D (circles), and 
dashed lines are guides to the eye. For comparison, we 
also present as solid lines the dependence of the vibra- 
tional energy of a single particle in a harmonic approxi- 
mation with the frequencies cjii and uj± given above. At 
low T the actual zero-point vibrational energy results 
to be smaller than that given in the harmonic approach 
[i.e., h(ui\\ + u>±/2)], but this trend changes as temper- 
ature rises. In fact, at temperatures larger than 1000 
K the vibrational energy derived from the simulations is 
larger than that obtained for the one-particle harmonic 
approach. 

Another way of getting insight into the anharmonic- 
ity of the defect vibrations is by looking at the ratio 
AE^/AE^ between vibrational energies of point defects 
for both isotopes. At low temperature, this ratio should 
converge to \[2 in the harmonic approximation, going 
to unity at high temperatures. The results yielded by 
our simulations are shown in Fig. 3, along with the one- 
particle harmonic expectancy. The former lie somewhat 
lower than the latter in the whole temperature region 
under consideration. In particular, at 300 K we find 
AEf/AE® = 1.364, vs a ratio of 1.394 derived in a har- 
monic approximation. 

PIMD simulations allow us to separate the potential 
{E p ) and kinetic (Ek) contributions to the vibrational 
energy. 42 ' 43 i 44 In fact, Ek is related to the spatial delo- 
calization of the quantum paths, which can be obtained 
directly from the simulations (see below). In Fig. 4, we 
present the kinetic and potential energy of the point de- 
fect associated to hydrogen, as a function of temper- 



ature. Symbols indicate results of our simulations for 
Ek (squares) and E p (circles), whereas dashed lines are 
guides to the eye. The potential energy of the point de- 
fect is found to be clearly larger than the kinetic energy, 
indicating an appreciable anharmonicity of the whole de- 
fect (in a harmonic approach one has Ek = E p ). A solid 
line represents the expected dependence for both, kinetic 
and potential contributions, in a harmonic approxima- 
tion with the frequencies uu and lu± given above. At 
low temperature we find an appreciable change of the 
kinetic energy respect to the harmonic approach, con- 
trary to the potential energy, which coincides within error 
bars with the harmonic expectancy. A qualitative under- 
standing of this behavior can be obtained by analyzing 
the energy changes obtained through time-independent 
perturbation methods. 22 ! Thus, assuming a perturbed 
one-dimensional harmonic oscillator (with perturbations 
of x 3 and x A type) at T — 0, the first-order change in 
the energy is totally due to a change in the kinetic en- 
ergy, the potential energy remaining unaltered respect to 
its unperturbed value. A similar behavior has been ob- 
tained for hydrogen in silicon from path-integral Monte 
Carlo simulations.— The main difference is that in that 
case the kinetic energy was found to increase respect to 
the harmonic value, contrary to the result obtained here 
for H on graphene. This seems to depend on the details 
of the interatomic interactions and the actual geometry 
of the point defect under consideration, but in both cases 
the potential energy at low temperature is very close to 
the value yielded by the harmonic approximation. 

Path-integral simulations at finite temperatures de- 
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FIG. 3: Temperature dependence of the ratio AE^/AE^ be- 
tween vibrational energies of hydrogen and deuterium defects. 
Symbols indicate results of PIMD simulations. For compari- 
son, we also present the ratio expected in a one-particle har- 
monic approximation (continuous line). The dashed line is a 
guide to the eye. 
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FIG. 5: Mean-square displacement D 2 of the quantum paths 
for hydrogen in directions parallel (circles) and perpendicular 
(squares) to the graphene sheet. Dashed lines correspond to 
a harmonic approximation for vibrations with frequencies uv 
and uj±_ . Error bars of the simulations results are on the order 
of the symbol size. 
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FIG. 4: Temperature dependence of kinetic (squares) and 
potential (circles) contributions to the vibrational energy of 
the H defect. A solid line represents the expectancy of a 
harmonic approximation with frequencies u)u andtJx- Dashed 
lines are guides to the eye. 



scribe quantum delocalization through paths of finite 
size. This means that the average extension of the paths 
is a measure of the importance of quantum effects in a 
given problem. One can define a kind of quantum de- 
localization as the mean-square "radius of gyration" D 2 
of the ring polymers associated to the quantum particle 



under consideration.-!^ 144 This means: 



(6) 



where x is the position of the centroid of the paths de- 
fined in Eq. JT]), and (...) indicates a thermal average 
at temperature T. Note that the total spacial delocal- 
ization of a particle includes an additional term, taking 
into account displacements of the center of gravity of 
the paths. This term is the only one surviving at high 
temperatures, since in the classical limit each path col- 
lapses onto a single point. For our problem of hydrogen 
on graphene, with well-defined H vibrations along three 
orthogonal axes, one can define a quantum delocaliza- 
tion of hydrogen along each of these directions (D 2 and 

.Dj_). These derealizations are displayed in Fig. 5 for vi- 
brations parallel (circles) and perpendicular (squares) to 
the graphene plane, as derived from our PIMD simula- 
tions. D 2 decreases as the temperature is raised and the 
particle becomes more "classical". For comparison, we 
also show in Fig. 5 the mean-square displacement D 2 ex- 
pected for harmonic oscillators of frequencies u>\\ and lj±, 
which can be worked out analytically 17 ' 44 i 46 The delo- 
calization D 2 derived from the simulations follows closely 
the harmonic expectancy in the whole temperature range 
considered here. However, D\ is higher than its corre- 
sponding harmonic result at T < 400 K. 
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B. Hydrogen diffusion 

Hydrogen is expected to diffuse on the graphene sheet, 
breaking a C-H bond and forming a new one with a 
nearby C atom. To check the appearance of this kind 
of hydrogen jumps, we have carried out some classi- 
cal molecular dynamics simulations, which allow one to 
follow the adatom along its trajectory on the surface. 
As indicated above, this kind of classical limit is eas- 
ily achieved with our PIMD code, by setting the Trotter 
number L = 1 (one bead per atom). At temperatures 
lower than 1000 K, hydrogen jumps on the graphene sur- 
face result to be infrequent events, and are rarely ob- 
served along a simulation run. Then, a reliable esti- 
mation of the diffusion coefficient by this method is not 
possible. Moreover, at temperatures larger than 1000 K 
hydrogen begins to escape from the sheet, without re- 
maining on it enough time for studying quantitatively 
the diffusion process. This situation is remedied by an- 
alyzing the hydrogen motion by TST, and in particular 
by the quantum version presented in Sect. II. B. This pro- 
cedure allows us to obtain free energy barriers for impu- 
rity jumping, which include corrections due to the quan- 
tum character of the atomic nuclei, and in particular the 
renormalization of the barriers caused by zero-point mo- 
tion. Also, phonon-assisted tunneling is included in the 
calculation, since motion of the C atoms takes into ac- 
count a full quantization of the vibrational degrees of 
freedom. The main question we address here is the de- 
pendence of hopping rates on temperature and impurity 
mass. In connection with this, there is a vast literature 
about theoretical models for quantum diffusion of light 
particles in solids. 17 ' 47 i 48 ' 49 Due to the complexity of this 
problem, such calculations have been typically based on 
model potentials for the impurity-lattice interactions. 

Here, the jump rate of the impurity between two near- 
est equilibrium sites is derived from the free-energy bar- 
rier between those sites. To calculate this barrier, we 
first select a continuous path from one site to the other, 
minimizing the energy at the transition (saddle) point in 
a classical calculation (point-like atoms) at T = 0. In 
connection with this, it is worthwhile noting that seem- 
ingly simple atomic jumps can actually involve coupled 
barriers, as indicated in Ref.[5(J Thus, to obtain the bar- 
rier for hydrogen jumps, we considered a coupled motion 
of H and the nearest C atom. For the optimal path, we 
obtained an energy barrier of 0.78 eV, with the saddle 
point corresponding to a symmetric (bridge) configura- 
tion of hydrogen between two C atoms and at a distance 
of 1.7 A from the graphene plane. Once selected the 
best path, we performed finite-temperature path-integral 
simulations with the centroid of H fixed on several points 
along this path, as described in Sect. II. B. Then, the free- 
energy barrier is calculated from the mean force by using 
Eq. ©. 

In Fig. 6 we present the free-energy barrier AF for 
adatom diffusion between adjacent C-H bonds. Data de- 
rived from line integration of the mean force are shown 
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FIG. 6: Effective free-energy barrier for impurity jumps on 
the graphene sheet as a function of temperature. Squares, hy- 
drogen; circles, deuterium; triangles, classical limit. Dashed 
lines are guides to the eye. 



as a function of temperature for hydrogen (squares) and 
deuterium (circles). For comparison, we also show re- 
sults derived from classical simulations (triangles). In 
this plot, one notices first that AF is higher for D than 
for H, and even higher for the classical limit, which in 
this respect may be considered as the large-mass limit 
M — ► oo i 33 i 35 Then, we observe an appreciable decrease 
in the free-energy barrier as the impurity mass is reduced. 
Second, one observes an increase in AF as temperature 
is raised in the three cases shown in Fig. 6. This increase 
is smaller, but not negligible, in the classical limit, which 
at T — > converges to the potential energy barrier given 
above (AE p = 0.78 eV). At low temperature, the depen- 
dence of AF upon impurity mass is related to the change 
in internal energy E of the defect complex along the dif- 
fusion path, since then the entropy contribution to the 
free energy becomes negligible. At 300 K we find AF = 
0.79, 0.74, and 0.71 eV, for the classical limit, D, and H, 
respectively. This means that quantum effects renormal- 
ize the free-energy barrier at room temperature by about 
6 % for D and 10 % for H. This reduction in AF decreases 
as temperature rises and the atoms become "more classi- 
cal". In fact, these free-energy barriers should converge 
one to the other in the high-temperature limit, which 
cannot be approached here due to the onset of adatom 
desorption. 

Shown in Fig. 7 is the rate for impurity jumps on a 
graphene sheet as a function of the inverse temperature. 
Data were derived from the free-energy barriers displayed 
in Fig. 6, by using Eq. ([2]). Results are given for hydrogen 
(squares), deuterium (circles), and classical limit (trian- 
gles). At T < 500 K, the jump rate for hydrogen results 
to be larger than that for deuterium, which in turn is 
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FIG. 7: Rate for impurity jumps on the graphene sheet as 
a function of the inverse temperature. Symbols represent re- 
sults derived from simulations for hydrogen (squares), deu- 
terium (circles), and classical limit (triangles). Dashed lines 
are guides to the eye. 



higher than that found in the classical limit, as expected 
from the change in effective free-energy barrier discussed 
above. At 300 K, the jump rate for hydrogen is found to 
be 11.6 s _1 , about 20 times larger than the value found 
in the classical calculation. The influence of quantum ef- 
fects on hydrogen diffusivity increases as temperature is 
lowered, as could be expected, but the jump rate itself 
becomes very small. In fact, at 200 K the calculated rate 
for hydrogen, /ch, is less than 10 -4 s -1 . On the contrary, 
at high temperatures fcn converges to the classical result, 
and the difference between both becomes unobservable 
at T larger than 1000 K. At this temperature we find 
fen = 2.0 x 10 9 s" 1 . 

The diffusion barriers obtained here for H and D on 
graphene are comparable to those found for hydrogen dif- 
fusion on graphite. Ferro et al£^ used density-functional 
theory to study various diffusion barriers in graphite. For 
hydrogen diffusion on its surface, they obtained a (classi- 
cal) barrier of 0.94 eV, somewhat larger than that found 
here at T = in the classical limit. 

The diffusion of H and D on solid surfaces has been 
studied earlier using quantum TST, employing simu- 
lations similar to those presented here. Such studies 
were mainly focused on hydrogen diffusion on metal 
surfaces ! 20 ' 52 Our results are qualitatively similar to 
those obtained in these studies, with an appreciable en- 
hancement of the jump rate of the impurity in com- 
parison to a classical model. At temperatures lower 
than those studied here (T < 50 K), a temperature- 
independent diffusion was found on metal surfaces ! 20 ' 53 
This is also a possibility for hydrogen diffusion on 
graphene, and remains as a challenge for future research. 
PIMD simulations at those temperatures would need the 



use of Trotter numbers L larger than 100, that together 
with the interaction potentials employed here requires 
computational resources out of the scope of the present 
work. 



IV. SUMMARY 

The main advantage of PIMD simulations of hydro- 
gen on graphene is the possibility of calculating defect 
energies at finite temperatures, including a full quantiza- 
tion of host-atom motions, which are not easy to take 
into account from fixed-lattice calculations and classi- 
cal simulations. Isotope effects can be readily explored, 
since the impurity mass appears as an input parame- 
ter in the calculations. This includes the consideration 
of zero-point motion, which together with anharmonic- 
ity may cause appreciable non-trivial effects. Our results 
indicate that hydrogen adsorbed on graphene cannot be 
accurately described as a particle moving in a harmonic 
potential. Even if anharmonicities of the interatomic po- 
tential are taken into account, a single-particle approx- 
imation is questionable as a realistic description of im- 
purity complexes at finite temperatures. It is then nec- 
essary to treat the defect as a many-body problem with 
anharmonic interactions. 

Ab-initio theoretical techniques to calculate defect en- 
ergies in solids have achieved an excellent precision in 
recent years. However, zero-point motion is a factor lim- 
iting the accuracy of state-of-the-art techniques to pre- 
dict energy bands and total energies of solids^ 5 - The same 
happens for defect levels caused by impurities in solids, 
since their energy may change appreciably as the impu- 
rity mass is varied. Here we have illustrated how an- 
harmonicities in the atomic motion cause an appreciable 
difference between kinetic and potential energy of the de- 
fect (about 20 % for H at 300 K), and have quantified 
the effect of the impurity mass on anharmonic shifts in 
the energy. 

Due to the large relaxation of the nearest C atoms, 
hydrogen migration requires important motion of these 
atoms. Then, an adatom jump has to be viewed as a 
cooperative process involving a coupled motion of the 
impurity and the nearest host atoms. This picture is 
similar to that described in the literature as "opening of a 
door" , 54 which favors impurity diffusion. Thus, quantum 
motion of both adatom and C atoms helps to rcnormalizc 
the diffusion barriers respect to the classical expectancy. 
At 300 K we find a hydrogen diffusivity 20 times larger 
than that derived from classical barriers. 



Acknowledgments 

This work was supported by Ministerio de Ciencia e 
Innovation (Spain) through Grant No. BFM2003-03372- 
C03-03. 



8 



1 A. K. Geim and K. S. Novoselov, Nat. Mater. 6, 183 (2007). 

2 M. I. Katsnelson, Mater. Today 10, 20 (2007). 

3 A. C. Dillon and M. J. Heben, Appl. Phys. A 72, 133 
(2001). 

4 M. H. F. Sluiter and Y. Kawazoe, Phys. Rev. B 68, 085410 
(2003). 

5 S. J. Pearton, J. W. Corbett, and M. Stavola, Hydrogen in 
Crystalline Semiconductors (Springer, Berlin, 1992). 

6 S. K. Estreicher, Mater. Sci. Eng. R14, 319 (1995). 

7 J. C. Meyer, C. O. Girit, M. F. Crommie, and A. Zettl, 
Nature 454, 319 (2008). 

8 R. Zeisel, C. E. Nebel, and M. Stutzmann, Appl. Phys. 
Lett. 74, 1875 (1999). 

9 O. V. Yazyev and L. Helm, Phys. Rev. B 75, 125408 
(2007). 

10 S. Casolo, O. M. Lovvik, R. Martinazzo, and G. F. Tanta- 
rdini, J. Chem. Phys. 130, 054704 (2009). 

11 D. W. Boukhvalov, M. I. Katsnelson, and A. I. Lichten- 
stein, Phys. Rev. B 77, 035427 (2008). 

12 E. J. Duplock, M. Schemer, and P. J. D. Lindan, Phys. 
Rev. Lett. 92, 225502 (2004). 

13 P. L. de Andres and J. A. Verges, Appl. Phys. Lett. 93, 
171915 (2008). 

14 F. Buda, G. L. Chiarotti, R. Car, and M. Parrinello, Phys. 
Rev. B 44, 5908 (1991). 

15 G. Panzarini and L. Colombo, Phys. Rev. Lett. 73, 1636 

(1994) . 

16 S. Bedard and L. J. Lewis, Phys. Rev. B 61, 9895 (2000). 

17 M. J. Gillan, Phil. Mag. A 58, 257 (1988). 

18 D. M. Ceperley, Rev. Mod. Phys. 67, 279 (1995). 

19 D. E. Makarov and M. Topaler, Phys. Rev. E 52, 178 

(1995) . 

20 T. R. Mattsson and G. Wahnstrom, Phys. Rev. B 51, 1885 
(1995). 

21 R. Ramirez and C. P. Herrero, Phys. Rev. Lett. 73, 126 

(1994) . 

22 C. P. Herrero and R. Ramirez, Phys. Rev. B 51, 16761 

(1995) . 

23 T. Miyake, T. Ogitsu, and S. Tsuneyuki, Phys. Rev. Lett. 
81, 1873 (1998). 

24 C. P. Herrero and R. Ramirez, Phys. Rev. Lett. 99, 205504 
(2007). 

25 M. C. Gordillo, J. Boronat, and J. Casulleras, Phys. Rev. 
B 65, 014503 (2001). 

26 M. C. Gordillo, Phys. Rev. B 76, 115402 (2007). 

27 R. P. Feynman, Statistical Mechanics (Addison- Wesley, 
New York, 1972). 

28 H. Kleinert, Path Integrals in Quantum Mechanics, Statis- 
tics and Polymer Physics (World Scientific, Singapore, 
1990). 

29 G. J. Martyna, M. E. Tuckerman, D. J. Tobias, and M. L. 
Klein, Mol. Phys. 87, 1117 (1996). 



M. E. Tuckerman, in Quantum Simulations of Complex 
Many-Body Systems: From Theory to Algorithms, edited 
by J. Grotendorst, D. Marx, and A. Muramatsu (NIC, FZ 
Jiilich, 2002), p. 269. 

31 D. Porezag, T. Frauenheim, T. Kohler, G. Seifert, and 
R. Kaschner, Phys. Rev. B 51, 12947 (1995). 

32 C. M. Goringe, D. R. Bowler, and E. Hernandez, Rep. 
Prog. Phys. 60, 1447 (1997). 

33 C. P. Herrero, R. Ramirez, and E. R. Hernandez, Phys. 
Rev. B 73, 245211 (2006). 

34 M. E. Tuckerman and A. Hughes, in Classical and Quan- 
tum Dynamics in Condensed Phase Simulations, edited by 
B. J. Berne, G. Ciccotti, and D. F. Coker (Word Scientific, 
Singapore, 1998), p. 311. 

35 R. Ramirez, C. P. Herrero, and E. R. Hernandez, Phys. 
Rev. B 73, 245202 (2006). 

36 G. A. Voth, D. Chandler, and W. H. Miller, J. Chem. Phys. 
91, 7749 (1989). 

37 J. Cao and G. A. Voth, J. Chem. Phys. 100, 5106 (1994). 

38 M. J. Gillan, J. Phys. C: Solid State Phys. 20, 3621 (1987). 

39 G. A. Voth, J. Chem. Phys. 97, 8365 (1993). 

40 C. P. Herrero, Phys. Rev. B 55, 9235 (1997). 

41 J. C. Noya, C. P. Herrero, and R. Ramirez, Phys. Rev. 
Lett. 79, 111 (1997). 

42 M. F. Herman, E. J. Bruskin, and B. J. Berne, J. Chem. 
Phys. 76, 5150 (1982). 

43 A. Giansanti and G. Jacucci, J. Chem. Phys. 89, 7454 
(1988). 

44 M. J. Gillan, in Computer Modelling of Fluids, Polymers 
and Solids, edited by C. R. A. Catlow, S. C. Parker, and 
M. P. Allen (Kluwer, Dordrecht, 1990), p. 155. 

45 L. D. Landau and E. M. Lifshitz, Quantum Mechanics 
(Pergamon, Oxford, 1965), 2nd ed. 

46 R. Ramirez and C. P. Herrero, Phys. Rev. B 48, 14659 

(1993) . 

47 C. P. Flynn and A. M. Stoneham, Phys. Rev. B 1, 3966 
(1970). 

48 H. Sugimoto and Y. Fukai, Phys. Rev. B 22, 670 (1980). 

49 H. R. Schober and A. M. Stoneham, Phys. Rev. Lett. 60, 
2307 (1988). 

50 M. Ramamoorthy and S. T. Pantelides, Phys. Rev. Lett. 
76, 267 (1996). 

51 Y. Fcrro, F. Martinelli, A. Jelea, and A. Allouche, J. Chem. 
Phys. 120, 11882 (2004). 

52 S. W. Rick, D. L. Lynch, and J. D. Doll, J. Chem. Phys. 
99, 8183 (1993). 

53 T. R. Mattsson and G. Wahnstrom, Phys. Rev. B 56, 
14944 (1997). 

54 D. E. Boucher and G. G. DeLeo, Phys. Rev. B 50, 5247 

(1994) . 



